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Abstract 

A double-atom partitioning of the molecular one-electron density matrix is used to describe atoms and 
bonds. All calculations are performed in Hilbert space. The concept of atomic weight functions (familiar 
from Hirshfeld analysis of the electron density) is extended to atomic weight matrices. These are con- 
structed to be orthogonal projection operators on atomic subspaces, which has significant advantages in 
the interpretation of the bond contributions. In close analogy to the iterative Hirshfeld procedure, self- 
consistency is built in at the level of atomic charges and occupancies. The method is applied to a test set 
of about 67 molecules, representing various types of chemical binding. A close correlation is observed 
between the atomic charges and the Hirshfeld-I atomic charges. 



1 



I. INTRODUCTION 

The concept of Atoms In the Molecule (AIM) has always played a central role in classifying 
and predicting chemical properties. Because this concept does not naturally show up in molecu- 
lar orbital (MO) theory, there is a sustained interest in extracting chemical atoms and functional 
groups from MO-based calculations. Among the most popular techniques are MuUiken popula- 
tion analysis yj, Bader's Quantum Chemical Topology (QCT) [12]-H1, the Hirshfeld method in its 
original ^ and iterative version [i6i-|9|, the iterative stockholder approach [ITOl . and Mayer's fuzzy 
atoms [fTTTl . 

Most methods are restricted to the partitioning of the electron density, but not all properties 
of a quantum mechanical object can be explicitly expressed in terms of the electron density. A 
more fundamental approach to the AIM should be based on density matrices [ fT2| . Because of the 
inherent non-locality of the density matrix, it was recently argued that a two-index partitioning 
into atomic (diagonal) and bond (off-diagonal) contributions, is necessary to guarantee the local 
nature and transferability of the atoms [[T3l[T4l . The authors recently introduced such an approach 
lfT4ll . but the bond matrices had significant contributions from core electrons and free electron 
pairs, which somewhat blurred the interpretation in terms of chemical bonding. In this paper we 
attempt to improve the description of atoms and bonds in a molecule. The interpretive problems 
are overcome by defining atomic weight matrices to be orthogonal projection operators onto one- 
electron subspaces assigned to atoms. The proposed method bears a close resemblance to the 
iterative Hirshfeld procedure, as we introduce the requirement of consistency between weight 
matrices and atomic contributions. 

n. DOUBLE ATOM PARTITIONING SCHEME 

The spin-summed one-electron density matrix (IDM) for a singlet iV-electron molecular wave- 
function ^{xi, . . . ,xiy) is expressed as 




(1) 



a 




dx2... / dxN'^\r'a,X2, ...XN)'^{ra,X2, ...xn) , 
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where the matrix representation in an orthogonal basis set (pi{r) was used. A decomposition of 
the identity matrix into atomic weight matrices [fTSl {wA)ij 

^ij = ^i^A)ij, (2) 

A 

can be inserted on both sides of the molecular one-electron density matrix to arrive at a double- 
index atomic partitioning 

P = ^ PAB = ^ WApWB- (3) 

AB AB 

The characteristics of such a partitioning are spatially localized contributions paa belonging to an 
atom A and some delocalized contributions {pab) b^a that describe the bonding between atoms A 
andB lfT4]| . 

If the atomic weight matrices are chosen as projectors onto orthogonal subspaces, which im- 
plies that they are both idempotent and orthogonal 

waWb = waSab, (4) 

then the bond matrices (pab) Bi^A have a zero trace and only the atomic matrices paa will con- 
tribute to the number of electrons N in the molecule 

N = Tr{p) = Y,TripAA). (5) 

A 

Note that in general the orthogonality of the weight matrices does not imply a zero overlap between 
the atomic electron densities PAAi^, r') in coordinate space, as occurs in e.g. QCT. 



ni. CONSTRUCTING ORTHOGONAL ATOMIC PROJECTORS 

The decomposition of one-electron space according to Eq. (|4]) can be done in myriad ways, so 
a way must be found that generates a chemically relevant decomposition. We found that, starting 
from the IDM's p^]^ for isolated atoms, the following recursive scheme {i = 0, 1, ...) 

P« = J:paa; ^ = {f^)-'^f^Af^)-^ 

A 

Paa = wypwy (6) 

converges and generates weight matrices Wa°^ that obey Eq. (4). In a forthcoming publication 
we will analyze this result and introduce alternative recursive schemes, including proofs. Here we 



just mention three points: (1) it is understood that the molecular IDM p used in Eq. ([6]) is positive 
semidefinite by construction and therefore the same holds for the and p^*^; (2) when the IDM 
of the isolated atom p^]^ contains only basis functions centered on atom A, the subspaces spanned 
by the eigenvectors of p^]^ are linearly independent, and as a result the weight matrices w'j^ of the 
0th iteration already obey Eq. Q; (3) depending on the level of theory for the atomic calculation it 
is possible that a nuUspace is generated for p^^\ making the inverse square root in Eq. ^ singular. 
We therefore replace any eigenvalues of the p^]^ and p by a small (10~^) positive value. This is 
sufficient to ensure that all p*-*^ remain nonsingular. 



rV. SELF-CONSISTENCY REQUIREMENTS 

As in ordinary Hirshfeld, the choice of the starting point (the isolated atom IDM's p^]^) is 
crucial, since different results are obtained depending on the charge and electronic state of the 
isolated atom. Following the ideas behind the iterative Hirshfeld procedure [61, the result can be 
made independent of the starting point by building in self-consistency through an outer iterative 
scheme (where Eq. (|6]) would represent the inner iterative scheme). Note that the need for a fitted 
start point 

In our simplest implementation (called 'charge equalization') we start from rotationally aver- 
aged IDM's of the neutral isolated atoms. The recursive scheme in Eq. ([6]) generates effective 
electron numbers Na = Tr(p^^) for the atoms. These can be used to create a rotationally av- 
eraged IDM of the charged isolated atom according to the linear interpolation between integer 
electron numbers (k < Na < A; + 1) 

P^IIINa] = {k + l-N^)pfA[k] 

+ iN^-k)pfAk + l]. (7) 

The charged PaaI-^a] be used as the next starting point in Eq. (6) and the whole process is 
repeated until convergence for the effective electron numbers. 

We noticed that, in contrast to the electron density, the IDM is much more sensitive to a mis- 
match in the orientation of one-electron orbitals in the molecule and in the isolated atom used to 
set up the AIM. In some cases, this resulted in rather large AIM charges. In order to solve this 
problem we also implemented a more sophisticated scheme (called 'population equalization'). 
The simple implementation of the previous paragraph is followed up to the point where p^^\ [^a] 
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is obtained. We then use the eigenvalue decompositions 

(oo) \ ^ (oo) (oo) (oo) /o\ 

Paa = 2^n\A,kVAA,kVAA,k (8) 

k 

PaaI^a] - 2^n^^iip^^i(p^^i (9) 

I 

to generate the new starting point to Eq. (|6]). This is given by Eq. (|9]), but with the rotationally 
averaged occupation numbers i replaced with the occupation numbers n^^j^ of the partitioned 
atom. The correspondence is made on the basis of maximal orbital overlap \{fAAk\fAAi)\- 
this procedure, the rotationally averaged IDM's are replaced by IDM's containing information on 
the molecular geometry. The whole process is again repeated until convergence. 



V. RESULTS AND DISCUSSION 



The procedure described in Sec. IV was tested by partitioning the IDM of a small set of ca. 67 
simple molecules with a singlet ground state, representative of a variety of chemical bonds. The set 
of molecules comprises the species tested in [14], supplemented with some extra molecules with 
relatively high Hirshfeld-I charges like CF4 and H2SO4. The molecular and atomic IDM's were 
calculated at the restricted Hartree-Fock level of theory using the Aug-cc-pVDZ basis set [[T6] - fT8l . 
The molecular geometry was taken from a B3LYP [I20ti23ll /cc-pVDZ [[nl [HI EH optimization. 
For CO additional calculations were performed in larger (Aug-cc-pVTZ and Aug-cc-pVQZ) 
[T8l basis sets in order to assess basis set convergence. 



A. The CO molecule 

Diagonalization of the matrices in Eq. ([3]) leads to natural orbitals and occupancies. Fig. [T] 
presents the dominant natural orbitals, and the corresponding occupancies, of the atomic density 
matrix of carbon (a-e) and oxygen (f-j) in a CO molecule using the population equalization scheme. 

Similar to the double-index partitioning in 3D space based on the use of Hirshfeld-I weights 
lfT4ll . the natural orbitals are slightly deformed versions of the typical atomic Is (a), 2s (b), 2pz, 2px 
and 2py (c-d-e) orbitals. However, in contrast to [|I4| . the current scheme generates occupancies 
of the orbitals not involved in bonding (core orbitals or orbitals containing free electron pairs) that 
are very close to their expected integer value. Also note that the double occupied molecular a -, tti 
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FIG. 1: The dominant natural orbitals, and the corresponding occupancies of the atomic density matrix of 
carbon (a-e) and oxygen (f-j) in a CO molecule. For details, see text. 
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FIG. 2: The dominant eigenvectors and eigenvalues of (twice) the CO bond matrix in a CO molecule. For 

details, see text. 



and -K2 Hartree-Fock orbitals simply divide their occupancy over the atomic pz, Px and Py orbitals. 
The summed occupancy of (c) and (h) in Fig. [T]is exactly 2. 

Fig. |2] depicts the dominant eigenvectors and their eigenvalues of the CO bond matrix pab in 
the CO molecule. Their shapes match those of the expected bonding and antibonding orbitals, 
but in contrast to [|14||. the corresponding bonding and antibonding orbitals have exactly opposite 
eigenvalues by construction. The a* a, ni* tti and tt2* ^2 electron relocations are the 
main contributions to the bond matrix, though there are some smaller contributions (i — )• d , j — )■ e) 
that have a mainly non-bonding character. Note that we avoid the term "occupancy" for the bond 
matrices, as the bond matrix reflects a change in density rather then a density itself. 

Apart from the dominant contributions shown in Figs.[T]j2} various orbitals with much smaller 
eigenvalues are also present (as the Aug-cc-pVDZ molecular basis set has 46 basis functions). 
A complete overview is given in Table |I] for both population equalization and charge equalization 
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C ,C 


0,0 


c,o + o,c 


(a) 


2.000 ( 2.000 ) 


2.000 ( 2.000 ) 


± 0.937 ( ± 0.919 ) 


(b) 


1.999 ( 1.988 ) 


1.987 ( 1.992) 


± 0.788 ( ± 0.799 ) 


(c) 


0.649 ( 0.605 ) 


1.615 ( 1.601 ) 


± 0.788 ( ± 0.799 ) 


(d) 


0.385 ( 0.399 ) 


1.615 ( 1.601 ) 


±0.163 (±0.154) 


(e) 


0.385 ( 0.399 ) 


1.350 ( 1.395 ) 


±0.037 (±0.128) 


(f) 


0.013 ( 0.008 ) 


0.001 ( 0.012 ) 


± 0.014 ( ± 0.009 ) 


(g) 


0.000 ( 0.000 ) 


0.000 ( 0.000 ) 


±0.011 (±0.007) 


(...) 


< 10-13 


< 10-13 


-10-13 X 10-13 



Sum 5.432 ( 5.400 ) 8.568 ( 8.600 ) 2.738 ( 2.816 ) 



TABLE I: All natural populations in the atomic density matrices (CC and 00) and eigenvalues of the bond 
matrix (CO) in a CO molecule. Non-bracketed and bracketed values belong to the population 
equaUzation scheme and the charge equalization scheme respectively. For details, see text. 

schemes (bracketed values). Note that for the atomic density matrices, apart from the five dominant 
natural orbitals, only two more have small (~ IQ-^ — 10-3) populations, while the remainder have 
vanishing (< 10 i3) populations. For the bond matrix, apart from the ten dominant orbitals, only 
four more have small eigenvalues, as was already noted in [fT4|. 

In Table [ll] the stability of the proposed density matrix partitioning is examined. The results 
seem to be quite stable with respect to basis set size, with differences for the traces of atomic 
density matrices less than 0.01 going from DZ to TZ, and less than 0.001 going from TZ to QZ. 
For the positive and negative components of the bond matrices there are fluctuations of about 
0.02 from DZ to TZ and from TZ to QZ. We checked that the stability holds even at the level of 
the individual orbital eigenvalues in Table |Ij In the charge equalization scheme, all eigenvalues 
differ less than 0.003 going from DZ to TZ, and less than 0.002 from TZ to QZ. In the population 
equalization scheme, the atomic matrices behave similarly (maximal deviation of 0.005 from DZ to 
TZ and of 0.003 from TZ to QZ), but the bond matrix eigenvalues deviate slightly more (maximal 
deviation of 0.015 from DZ to TZ and of 0.022 from TZ to QZ). We also checked that for almost 
all the orbitals in Figs. Figs.[T]j2} the shapes are visually indistinguishable when changing basis 
set size from DZ to QZ. The exceptions are the bond matrix orbitals (e) and (j), corresponding to 
small (in absolute value) eigenvalues, where a somewhat more distorted shape is observed. 
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CO 


Aug-cc-pVDZ 


Aug-cc-pVTZ 


Aug-cc-pVQZ 


A,B 


Tr(pA,B) 


Tr(PA,B) 


Tr(/9A,B) 


C,C 
0,0 

c,o + o,c 


5.432 ( 5.400) 
8.568 ( 8.600) 
± 2.738 ( ±2.816) 


5.425 ( 5.393) 
8.575 ( 8.607) 
± 2.756 (±2.811) 


5.426 (5.394) 
8.574 (8.606) 
±2.731 (±2.814) 



TABLE II: Basis set convergence of the summed positive and negative eigenvalues in the atomic and bond 
matrices for CO. Non-bracketed and bracketed values belong to the population equalization 
scheme and the charge equalization scheme respectively. For details, see text. 

B. Evaluation of the atomic charges 

Fig. [3] displays the correlation between the charges obtained with the Hirshfeld-I method ^ 
(known to be a reliable AIM method) and the charges resulting from either the charge equaliza- 
tion scheme or the population equalization scheme applied to the entire set of 67 molecules. Both 
schemes show a strong linear correlation with Hirshfeld-I. It is remarkable that the linear correla- 
tion is more satisfying for the population equalization scheme (R^=0.96, slope=1.18) than when 
only the atomic charges are made self-consistent (R^=0.90, slope = 1.45). As mentioned, this is 
related to the inadequacy of the rotationally averaged IDM's used as starting point in the charge 
equalization scheme. In some cases, the latter leads to unexpectedly high charges. This problem 
seems to be solved using the population equalization procedure, which is more adapted to the 
molecular geometry. 

VI. CONCLUSIONS 

We have implemented a two-index partitioning of the molecular density matrix into atomic and 
bond contributions using atomic weight matrices that are orthogonal projection operators in one- 
electron Hilbert space. The method is highly efficient in terms of computation time since no 3D 
integrals need to be computed. The resulting decomposition provides a rather physical descrip- 
tion of the adaptations made by atoms in a molecule and the deformation caused by the bonding 
process. The bond matrices are traceless, i.e. the electrons are all in the atomic contributions and 
not in the bond matrices. Core orbitals and free electron pairs are fully assigned to the atoms. The 
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New charges compared to Hirsh1eld-I charges 

4 I , , , 1 , 1 




-1 1 

Hirshfeld-I charges 



FIG. 3: Comparison of the Hirshfeld-I charges with those generated by the charge equalization scheme 
(red circles) and the population equalization scheme (black squares). For details, see text. 



trial atoms and the AIM atoms are required to have equal charges or equal orbital populations. 
The latter approach leads to a better correlation of the atomic charges with the Hirshfeld-I atomic 
charges. 
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